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We perform a lattice QCD calculation of the hadronic light-by-light scattering amplitude in a broad 
kinematical range. At forward kinematics, the results are compared to a phenomenological analysis 
based on dispersive sum rules for light-by-light scattering. The size of the pion pole contribution 
is investigated for momenta of typical hadronic size. The presented numerical methods can be 
used to compute the hadronic light-by-light contribution to the anomalous magnetic moment of the 
muon. Our calculations are carried out in two-flavor QCD with the pion mass in the range of 270 
to 450 MeV, and contain so far only the diagrams with fully connected quark lines. 


I. INTRODUCTION 

Light-by-light scattering, the elastic scattering of two 
photons, is a striking prediction of Quantum Electrody¬ 
namics (QED). The light-by-light (LbL) interaction ap¬ 
pears prominently in corrections to the anomalous mag¬ 
netic moment (g — 2) of the electron and muon. The 
muon (g — 2) exhibits a 3cr discrepancy between experi¬ 
ment and the Standard Model calculations [1]. While the 
current theory and experimental errors are comparable in 
size, a new (g — 2) M experiment [2] aiming to reduce the 
experimental error by a factor of four is in preparation 
at Fermilab. 

The theory error on (g — 2)^ is dominated by hadronic 
contributions, namely the hadronic vacuum polarization 
(HVP) and hadronic light-by-light (HLbL) scattering. 
Using unitarity and causality, the HVP contribution is 
expressed in terms of the total e + e~ —> hadrons cross 
section, and hence its precision can systematically be im¬ 
proved by collider experiments alone. By contrast, the 
HLbL contribution cannot be expressed entirely in terms 
of cross sections for 77-fusion into hadrons; see [3-5] for 
dispersive approaches to the problem. A direct ab initio 
calculation within Quantum Chromodynamics (QCD) is 
very challenging due to its non-perturbative nature. In 
this work we address the problem using lattice QCD. 

A first lattice QCD+QED calculation of the HLbL con¬ 
tribution to (g— 2) m has recently been performed by Blum 
et al. [6]. We envisage a different method where the four- 
point function for LbL scattering is computed in lattice 
QCD and integrated over to yield the HLbL contribution. 
In this Letter we present the four-point function calcu¬ 
lation and check it against the available phenomenology. 
Exploiting unitarity and causality, the forward HLbL am¬ 
plitude can be expressed as a dispersive integral over the 
7*7* —> hadrons cross section [7, 8]. A parametrization 
of the latter allows us to confront the lattice calculation 
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with phenomenology in a fairly straightforward manner. 
As the neutral pion (7T°) pole dominates the HLbL con¬ 
tribution to (g — 2)^ in phenomenological calculations [1], 
we study its relative size both at forward and off-forward 
kinematics. 

II. THEORY BACKGROUND 

The Lehmann-Symanzik-Zinnnermann reduction for¬ 
mula for the HLbL scattering amplitude implies [9] 

Mn llMli 2 fl 3 (pi,P 4 P 2 ,P 3 ) (1) 

= e 4 (—«n MlM2MAt4 (— p 4 ;-pi,p 2 )), 

where p^ = pi + P 4 — P 2 and 

in(p^PGP2) = j d 4 Xid 4 X 2 d' l X i (2) 

e +i Pa ' Xa (01T {j Ml (zi) j M2 O2) jp 3 (0)j M 4 (ar)} 10) 

is the Minkowski-space time-ordered correlator of the 
conserved vector current j^ = | xi'y^u— ^dj^d-i-.... The 
index a takes the values 1, 2 and 4. The components 
of the current used in the Euclidean theory [10] are 
related to their Minkowskian counterparts by Jo = jo, 
Jk = ijk- The analytic continuation then yields the fol¬ 
lowing relation to the Euclidean correlation function, 

-iIW 2 „ sM4 (HP?, P 4 ); HA 0 , Pi), (~iPl P2)) ( 3 ) 

= *" on m™4( p 4;Pi,P 2 ), 

n ™„4 ( p 4! Pi, Pi) = J d 4 x, d 4 X 2 d 4 X 4 (4) 

where no is the number of temporal indices carried by 
the vector currents in the correlator. 

The forward scattering case is obtained in Ecp (1) by 
setting p 2 = pi. Renaming the momenta to match the 
conventional notation, we have 

■^^2/43^4(91 >92) = A 4 MlP2AJ3M4 (gi,g 2 -A qi,q 2 ) ( 5 ) 

= e 4 (-'jIW 3 /M^ 2 (-9 2 ; -qi,qi))- 
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The forward scattering amplitude can be decomposed 
into eight Lorentz-invariant amplitudes [11]. They are 
functions of the virtualities qf and q 2 of the photons, as 
well as of the variable v = <7 ■ q 2 . Using the projector 
R^ onto the subspace orthogonal to qi and <; 2 , we focus 
here on the amplitude [12] 

MMvUlv) = (6) 


Combining Eqs. (5) and (3), we can access the amplitude 
A4 tt from the Euclidean correlator, 


A1tt(— Qi Q 2 ,— Qi ' Q2) 


R » v = ^" (<Qi ■ Q2) 2 - Q?Q 2 2 ' 

(Qi * Q2)( K QifiQ2u T Q11/Q2 74) 
-~QiQ2^Q-2v — QiQinQh 


( 7 ) 


( 8 ) 


The largest value of \u\ that can be reached with Eu¬ 
clidean kinematics is limited by the virtualities of the 
photons [13], \v\ < [QlQ 2 ) 1 ^ 2 < |(Qi +Ql) = Vo, while 
the nearest singularity is the s-channel 7r° pole located 
at tv = 2 {m 2 + Qf + Q 2 ). A technical issue arises 
when Qi and Q 2 are collinear: the projector R^ v be¬ 
comes ambiguous. To resolve the issue, we note that 
R[uj — R^ UifiUhi, where R^v = 6 ^ 1 , Qi^Qiiz/Qi 
and U\ is the unit vector parallel to the projection of 
Q 2 onto the subspace orthogonal to Q\. The average of 
the applied projector over the directions of U\ in that 
subspace yields 


U 1 — 5^711/4377/427*4 
Tyg ^7?7 ll yi 2 7?71 3 7l 4 + 


( 9 ) 


We use this averaged projector in Eq. (7) when Qi and 
Q 2 are collinear. 

In [8], it was shown that the HLbL amplitude A4tt{v), 
for fixed spacelike photon virtualities, can be obtained 
from the following dispersive sum rule, 


MTr(ql,ql,v) - A4 T t(<2i, q%, 0) 

2v2 f°° j < \fv przr q M 7 , u ^ 

= - / dv -TTT2 -2- + <J 2 ){v ), 

7T J„ 0 V'[V-V 1 - It) 


( 10 ) 
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FIG. 1. Four-point function quark contraction topologies. 
The vertices represent vector currents and the lines are quark 
propagators. In this work, we compute only the leftmost, 
fully-connected class of diagrams. 



FIG. 2. Fully-connected four-point function quark contrac¬ 
tions. Each panel represents two contractions with oppo¬ 
site directions of quark flow. The solid quark lines are com¬ 
puted using a point-source propagator, the dashed lines using 
sequential propagators, and the dotted lines using double- 
sequential propagators. 


with R(q 2 ,q 2 ) the pion transition form factor as defined 
in [14]. For q 2 = 0, the same result is obtained from the 
sum rule, using the expression for the 77* —> 7r° cross- 
section given in [8]. 

In summary, the amplitude A4 tt can be computed on 
the lattice via Eq. (7) and from e + e _ collider data via 
Eq. (10). In the following, we present a comparison of 
the two approaches. 

III. IMPLEMENTATION OF THE EUCLIDEAN 

FOUR-POINT FUNCTION IN LATTICE QCD 

In numerical lattice QCD calculations of n-point func¬ 
tions, the quark path integral is evaluated analytically to 
yield a sum of contractions of quark propagators. For the 
four-point function of vector currents, these fall into five 
distinct topologies, illustrated in Fig. 1. In this work, we 
compute only the six contractions that are fully quark- 
connected. 

We use a Wilson-type quark action, three lattice con¬ 
served currents J° and one site-local current J^ (see for 
instance [15] for an explicit definition). Generically, we 
evaluate the fully-connected contribution to 


where a 0 and cr 2 are the total cross sections 
7 * { ( li) 7 * {(l 2 ) hadrons with total helicity 0 and 2 re¬ 
spectively. It can be shown [8] that A4 tt vanishes at 
1 / = 0 if either of the photons is real. It is interesting to 
test the sum rule for the 7r° pole contribution. Using the 
expression for II Ml/p(T given in [14] and Eqs. (5, 6), one 
finds 


A4xt( — Q i> Q 2 R) 
?{-Q\-Q\Y 


= e 4 w 2 - QlQl) 
Q 1 +Q 2 + m 2 
{Q\ + Q\ + nn\) 2 - Av 2 


( 11 ) 


11 ^/ 427 * 3/44 (* 4 ; h,h)= E /l(*l)/ 2 (* 2 ) 
x lt x 2 

(j^ ( x i ( X 2 ) J ls (9) J/j ,4 (AA) + contact terms), (12) 

for some fixed functions /i. 2 and all values of {/x a } 
and X 4 . The contact terms are present when two or 
three lattice conserved currents coincide, and serve to 
ensure that the conserved-current relations hold, e.g., 
A{ t f 4) n|f i t M2/ 4 3 7 t4 = 0, where A X) is the backward lat¬ 
tice derivative. 
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x i () -5 m_ = 324 MeV, Of = 0.377 GeV 2 



Ql (GeV 2 ) 


FIG. 3. The forward scattering amplitude AItt at a fixed 
virtuality Qf = 0.377GeV 2 , as a function of the other photon 
virtuality Q \, for different values of v. The curves represent 
the predictions based on Eq. (10), see the text for details. 


The fully-connected contribution to Eq. (12) is evalu¬ 
ated using the method of sequential propagators. First, 
a point-source propagator is computed from X 3 . Then, 
it is combined with the function fi or /2 to form the 
source for a new (sequential) propagator. These sequen¬ 
tial propagators are then used to form sources for double- 
sequential propagators that depend on both fi and /2. 
Finally, the fully-connected contraction is formed using 
all three kinds of propagators; this is illustrated in Fig. 2. 
For generic complex fi and /2, this requires one point- 
source, 16 sequential and 32 double-sequential propaga¬ 
tors, although these counts can be reduced in various spe¬ 
cial cases. We have verified that in our implementation 
the four-point function matches the lattice perturbation 
theory calculation if the gauge link variables are set to 
unity, and that the conserved-current conditions hold on 
each gauge configuration. 

For evaluating the momentum-space correlator, we set 
the functions to be plane waves, / a (X) = e ~ lPa ' x and 
compute the Fourier modes with respect to X 4 . Thus, 
H^iC^i -Pi, P 2 ) can be evaluated efficiently at fixed 
Pi 2 for all P4 available on the lattice. 


IV. RESULTS 

We have used three lattice QCD ensembles with two 
degenerate flavors of non-perturbatively 0(a) improved 
Wilson quarks and a plaquette gauge action. The en¬ 
sembles are at a single lattice spacing a = 0.063fm [16], 
correspond to pion masses m n = 451, 324 and 277 MeV, 
and are respectively of spatial linear size 32, 48 and 48, 



FIG. 4. The dependence of the amplitude AItt on v, both 
photon virtualities being fixed at 0.377 GeV 2 , at three dif¬ 
ferent pion masses. The dashed and dotted curves show the 
7 T° and 7 r° + 7/ contributions (there is no 7 meson in two- 
flavor QCD), the solid curve includes all single-meson and 
7r + 7r _ contributions, and the dash-dotted curves additionally 
include the high-energy contribution for the case of real pho¬ 
tons at the physical pion mass. 


the time direction being twice as long; see [17] for more 
details. Only the up and down quark contributions to 
the electromagnetic current are included. The local vec¬ 
tor current is renormalized non-perturbatively [18]. 
The results shown here were obtained using fairly low 
statistics, with a maximum of 300 samples. 

Due to the finite volume of the lattice, the momenta 
take discrete values. The subtracted forward scatter¬ 
ing amplitude, M T t(~Qi, ~Qh ~Q\-, 0) 

(which is even in v), is obtained by linearly interpolating 
the second term between the available Q\ to match the 
first term. It is shown in Fig. 3 at fixed pion mass and 
fixed Qi, and also in Fig. 4 with both photon virtualities 
fixed. For the latter, linear interpolation in Q\ was also 
used in the first term, except for the points at maximal 
v. At fixed v, the amplitude tends to decrease as the 
virtualities are increased, at fixed virtualities it tends to 
increase with \v\, and at fixed kinematics we do not find 
a significant dependence on the pion mass. 

We compare the lattice data with results from the 
sum rule, Eq. (10), using a phenomenological model for 
the transverse 7*7* —► hadrons cross section, cto + cr 2 , 
based on Ref. [19]. We include pseudoscalar, scalar, 
axial-vector, and tensor mesons, as well as the non¬ 
resonant 7r + 7r _ contribution (in scalar treelevel QED 
with pion electromagnetic form factors). The 7*7* —> 
meson form factors have not been measured experi¬ 
mentally; they are assumed to factorize as F(qf,q 2 ) = 
F(ql,0)F(0, g|)/F(0, 0). For the pseudoscalar and axial- 
vector mesons, F(q 2 , 0) = F(0,q 2 ) is described based on 
experimental data as in Ref. [8] and, lacking guidance 
from experiment, we assume a monopole form factor for 
the scalar and tensor resonances with a pole mass set by 
hand to A = 1.6 GeV. The model is modified for unphys- 
































4 



FIG. 5. Two Lorentz contractions of the vector four-point 
function at non-forward kinematics. For A = I (squares), the 
pion pole contribution vanishes, while for A = 0 (circles), it 
does not. The curves correspond to the n° pole contribution 
in the latter case. 

ical quark masses by adjusting the masses and 77 decay 
widths, r 77 , of the mesons. The pion mass and decay 
constant [20] are calculated on each lattice ensemble, 
and .F(0,0) is set to the value (—inspired by 
the chiral anomaly prediction (see e.g. [21]). For each of 
the remaining mesons, the mass is assumed to have the 
same shift as that of the p meson, relative to the physi¬ 
cal point, and T 77 is assumed to scale linearly with the 
meson’s mass. 

This model together with the dispersive sum rule pro¬ 
duces the solid curves in Figs. 3 and 4, which agree well 
with the data. Varying A by ±0.4 GeV shifts the curves 
by up to ±50%, hence it is clear that the model has a con¬ 
siderable uncertainty; nevertheless the consistency with 
the data is remarkable. Fig. 4 also shows the individual 
contributions from 7r° and rf mesons and a high-energy 
contribution arising from a fit to the total 77 —> hadrons 
cross section [22] based on Regge theory. The latter is 
excluded from the main model curves due to the lack 
of a well-motivated extrapolation to the case of virtual 
photons and larger-than-physical pion masses. It is inter¬ 
esting to note that the two-pion production is typically 
the dominant contribution to the amplitude, rather than 
the 7r° and rf production. 

Moving to off-forward kinematics, the situation is more 
complicated. In general, the four-point function of vec¬ 
tor currents can be decomposed into 41 Lorentz-invariant 
functions [23] (see also [5]) that depend on six kinematic 
variables, of which three are fixed when Pi and P 2 are 
fixed in our lattice calculation. To study the importance 
of the 7r° contribution, we consider two contractions: 
n^ iMi/i3At3 , which has pion poles when (Pi ±P 4 ) 2 = -to 2 
or (P2 ±P 4 ) 2 = —to 2 , and a fully-symmetric contraction 
that has no 7r°-exchange contribution. These are shown 


in Fig. 5 , where we have also fixed Pf = P 2 [where 
P3 = —(Pi ± P2 ± P 4 )] to be a typical hadronic scale 
below 1 GeV 2 . We find that the fully-symmetric con¬ 
traction yields larger data, again indicating that the 7r° 
does not provide the dominant contribution. 

V. CONCLUSION 

We have demonstrated that the fully connected con¬ 
tribution to the momentum-space four-point function of 
the electromagnetic current can be computed with mo¬ 
derate computational effort in lattice QCD if two of the 
three momenta are fixed. As an application, we com¬ 
puted one of the forward 7*7* scattering amplitudes in 
a broad kinematic range. Via a dispersive sum rule, it 
is related model-independently to 7*7* —» hadrons cross 
sections. Modelling the latter, we find the comparison 
of the lattice calculation with the phenomenological ap¬ 
proach to be successful. The systematic uncertainties of 
the comparison are presently still large, mainly because 
our current calculations are performed at heavier quark 
masses than the physical ones, but this model depen¬ 
dence can be systematically reduced. Also, the not fully 
connected contraction topologies depicted in Fig. 1 could 
be important. We investigated the size of the pion pole 
contribution both in the forward and the off-forward am¬ 
plitude. Both the lattice data and the model show that 
it is by no means dominant in a range of kinematic in¬ 
variants of typical hadronic size. 

The numerical methods presented can be applied to 
a direct lattice calculation of the HLbL contribution to 
(g — 2) M : we are currently working on a position-space 
approach where the photon propagators are integrated 
out semi-analytically in infinite volume. The dominant 
systematic effects are likely to be quite different from 
those in the method of Blum et al. [6], allowing for useful 
cross-checks. Since phenomenological calculations indi¬ 
cate that the 7r° is dominant in the HLbL contribution 
to (g — 2)^ [1], realistically light quark masses and large 
volumes will be required to treat this long-range contri¬ 
bution correctly. Lattice data on the HLbL amplitude it¬ 
self can also help discriminate between phenomenological 
models used in the calculation of (g — 2) M . 
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